################################################################
# Replication File: "The Politics of Too Much - The Decline in Labor Power 
# and the Global Rise of Corporate Saving" 
# Author: Nils Redeker
# TSCS Analysis Replication
# R version 3.6.2 (2019-12-12)
################################################################

# clean working environment
rm(list = ls())


# load packages --- --- ---- --- --- --- --- --- --- --- --- --- -- --- --- --- --- --- --- -- --- --- --- --- --- 

library(ggplot2)
library(stargazer)
library(dplyr)
library(grid)
library(Amelia)
library(readr)
library(plm)
library(Hmisc)
library(ggpmisc)
library(pcse)
library(prais)
library(panelAR)
library(visreg)
library(cowplot)
library(mice)
library(rsq)
library(mitools)

# set working directory
setwd()

# load data --- --- ---- --- --- --- --- --- --- --- --- --- -- --- --- --- --- --- --- -- --- --- --- --- ----- 

df<-read_csv("~/Cross_nationa_data.csv") 


# --- --- ---- --- --- --- --- --- --- --- --- --- -- --- --- --- --- --- --- -- --- --- --- --- --- --- --- ---
# 1) FIGURE 1: DEVELOPMENT OF CORPORATE SAVING IN MAJOR ECONOMIES  --- --- ---- --- --- --- --- --- --- --- ----
# --- --- ---- --- --- --- --- --- --- --- --- --- -- --- --- --- --- --- --- -- --- --- --- --- --- --- --- ---

ggplot(df %>% filter(country %in% c("United States","United Kingdom","Germany","Italy", "France")), 
       aes(year, std_net_lending)) + geom_line(aes(linetype=country)) + geom_hline(yintercept = 0) +
  theme_minimal() + xlab("") + ylab("Corporate Net Lending (% of GDP)") +
  theme(legend.title = element_blank())


# --- --- ---- --- --- --- --- --- --- --- --- --- -- --- --- --- --- --- --- -- --- --- --- --- --- --- -- --- ---
# 2) DATA IMPUTATION WITH AMELIA  --- --- ---- --- --- --- --- --- --- --- --- --- --- --- -- --- --- --- --- --- -
# --- --- ---- --- --- --- --- --- --- --- --- --- -- --- --- --- --- --- --- -- --- --- --- --- --- --- -- --- --- 

df<-as.data.frame(df) # turn into data frame 

set.seed(12345) # set seed for replication of imputaition

a.out<-amelia(df, ts="year", cs="country", intercs = TRUE, m=5) # conduct 5 imputations 


#------------------------------------------------------------------------------------------------------------------------
# 3) TABLE 1 - TSCS: CORRELATION BETWEEN TRADE UNION DENSITY/ EMPLOYMENT PROTECTION AND CORPORATE SAVING
#------------------------------------------------------------------------------------------------------------------------

### 3.1) Regression Models Trade Union Density 

for(i in 1:a.out$m) {}

## Table 1 - Model 1

# Conduct regressions across all imputations 
mod1_ud <-lapply(a.out$imputations, function(x){
  panelAR(std_net_lending ~  ud + factor(country), 
          data = a.out$imputations[[i]],
          timeVar ="year",
          panelVar="country", autoCorr="psar1",
          panelCorrMethod="pcse")
})

# Extract standard errors
mod1_ud_se <- list()
mod1_ud_se <- lapply(1:5, function(x){
  sqrt(diag(mod1_ud[[x]]$vcov)) 
})

# Take averages of standard errors across the models 
mod1_ud_se <-as.numeric(Reduce("+", mod1_ud_se)/5)

# Extract coefficients
mod1_ud_coef <- list()
mod1_ud_coef <- lapply(1:5, function(x){
  mod1_ud[[x]]$coefficients
})

# Take averages of coefficients across the models 
mod1_ud_coef <- Reduce("+", mod1_ud_coef)/5

## Table 1 - Model 2

# Conduct regressions across all imputations 
mod2_ud <-lapply(a.out$imputations, function(x){
  panelAR(std_net_lending ~  ud + rti_score +  FDI_standardized + realgdpgr + 
            realinterest +  stock_market + old_age + cit + 
            + factor(country) + factor(year), 
          data = a.out$imputations[[i]],
          timeVar ="year",
          panelVar="country", autoCorr="psar1",
          panelCorrMethod="pcse")
})

# Extract standard errors
mod2_ud_se <- list()
mod2_ud_se <- lapply(1:5, function(x){
  sqrt(diag(mod2_ud[[x]]$vcov)) 
})

# Take averages of standard errors across the models 
mod2_ud_se <-as.numeric(Reduce("+", mod2_ud_se)/5)

# Extract coefficients
mod2_ud_coef <- list()
mod2_ud_coef <- lapply(1:5, function(x){
  mod2_ud[[x]]$coefficients
})

# Take averages of coefficients across the models 
mod2_ud_coef <- Reduce("+", mod2_ud_coef)/5

## Table 1 - Model 3

# Conduct regressions across all imputations 
mod3_ud <-lapply(a.out$imputations, function(x){
  pcse(lm(std_net_lending ~  ud  + rti_score + FDI_standardized  + realgdpgr + realinterest  +
            stock_market  + old_age + cit + std_net_lending_lag + ud_lag +  
            factor(country), data = a.out$imputations[[i]]),
       groupT =a.out$imputations[[i]]$year, groupN=a.out$imputations[[i]]$country)
})

# Extract standard errors
mod3_ud_se <- list()
mod3_ud_se <- lapply(1:5, function(x){
  mod3_ud[[x]]$pcse
})

# Take averages of standard errors across the models 
mod3_ud_se <-as.numeric(Reduce("+", mod3_ud_se)/5)

# Extract coefficients
mod3_ud_coef <- list()
mod3_ud_coef <- lapply(1:5, function(x){
  mod3_ud[[x]]$b
})

# Take averages of coefficients across the models 
mod3_ud_coef <- as.numeric( Reduce("+", mod3_ud_coef)/5)



### 3.2) Regression Models Employment Protection -----------------------------------------------------------------
for(i in 1:a.out$m) {}

## Table 1 - Model 4

# Conduct regressions across all imputations 
mod1_emprot_reg <-lapply(a.out$imputations, function(x){
  panelAR(std_net_lending ~  emprot_reg + factor(country), 
          data = a.out$imputations[[i]],
          timeVar ="year",
          panelVar="country", autoCorr="psar1",
          panelCorrMethod="pcse")
})

# Extract standard errors
mod1_emprot_reg_se <- list()
mod1_emprot_reg_se <- lapply(1:5, function(x){
  sqrt(diag(mod1_emprot_reg[[x]]$vcov)) 
})

# Take averages of standard errors across the models 
mod1_emprot_reg_se <-as.numeric(Reduce("+", mod1_emprot_reg_se)/5)

# Extract coefficients
mod1_emprot_reg_coef <- list()
mod1_emprot_reg_coef <- lapply(1:5, function(x){
  mod1_emprot_reg[[x]]$coefficients
})

# Take averages of coefficients across the models 
mod1_emprot_reg_coef <- Reduce("+", mod1_emprot_reg_coef)/5


## Table 1 - Model 5

# Conduct regressions across all imputations 
mod2_emprot_reg <-lapply(a.out$imputations, function(x){
  panelAR(std_net_lending ~  emprot_reg + rti_score +  FDI_standardized + realgdpgr + realinterest +  stock_market + old_age + cit + 
            + factor(country) + factor(year), 
          data = a.out$imputations[[i]],
          timeVar ="year",
          panelVar="country", autoCorr="psar1",
          panelCorrMethod="pcse")
})

# Extract standard errors
mod2_emprot_reg_se <- list()
mod2_emprot_reg_se <- lapply(1:5, function(x){
  sqrt(diag(mod2_emprot_reg[[x]]$vcov)) 
})

# Take averages of standard errors across the models 
mod2_emprot_reg_se <-as.numeric(Reduce("+", mod2_emprot_reg_se)/5)

# Extract coefficients
mod2_emprot_reg_coef <- list()
mod2_emprot_reg_coef <- lapply(1:5, function(x){
  mod2_emprot_reg[[x]]$coefficients
})

# Take averages of coefficients across the models 
mod2_emprot_reg_coef <- Reduce("+", mod2_emprot_reg_coef)/5


## Table 1 - Model 5
# Conduct regressions across all imputations 
mod3_emprot_reg <-lapply(a.out$imputations, function(x){
  pcse(lm(std_net_lending ~  emprot_reg  + rti_score + FDI_standardized  + realgdpgr + realinterest  +
            stock_market  + old_age + cit + std_net_lending_lag +  emprot_reg_lag + factor(country), data = a.out$imputations[[i]]),
       groupT =a.out$imputations[[i]]$year, groupN=a.out$imputations[[i]]$country)
})

# Extract standard errors
mod3_emprot_reg_se <- list()
mod3_emprot_reg_se <- lapply(1:5, function(x){
  mod3_emprot_reg[[x]]$pcse
})

# Take averages of standard errors across the models 
mod3_emprot_reg_se <-as.numeric(Reduce("+", mod3_emprot_reg_se)/5)

# Extract coefficients
mod3_emprot_reg_coef <- list()
mod3_emprot_reg_coef <- lapply(1:5, function(x){
  mod3_emprot_reg[[x]]$b
})

# Take averages of coefficients across the models 
mod3_emprot_reg_coef <- as.numeric( Reduce("+", mod3_emprot_reg_coef)/5)


### 3.3) Prepare "container" models for table 1 -----------------------------------------------------------------

# Note: these are placeholder models to insert coeffiecents/ standard errors estimates above into an easy table
# with stargazer

# Table 1 - Model 1 - placeholder
container_basis<- lm( std_net_lending ~ ud + factor(country), 
                      data = a.out$imputations$imp1)  
# Table 1 - Model 2 - placeholder
container_twoway<- lm(std_net_lending ~  ud + rti_score +  FDI_standardized + 
                        realgdpgr + realinterest +  stock_market + old_age + cit + 
                       + factor(country) + factor(year), data = a.out$imputations$imp1)
# Table 1 - Model 3 - placeholder
container_lag<-lm(  std_net_lending ~  ud  +  rti_score +   FDI_standardized +
                      realgdpgr +  realinterest  +
                      +  stock_market +  old_age +  cit +
                      +  std_net_lending_lag + ud_lag + factor( country), 
                    data = a.out$imputations$imp1) 
# Table 1 - Model 4 - placeholder
container_basis_emprot<- lm( std_net_lending ~ emprot_reg + factor(country), 
                      data = a.out$imputations$imp1)
# Table 1 - Model 5 - placeholder
container_twoway_emprot<- lm(std_net_lending ~  emprot_reg + rti_score +  FDI_standardized + 
                               realgdpgr + realinterest +  stock_market + old_age + cit + 
                        + factor(country) + factor(year), data = a.out$imputations$imp1)
# Table 1 - Model 6 - placeholder
container_lag_emprot<-lm(  std_net_lending ~  emprot_reg  +  rti_score +   FDI_standardized +
                      realgdpgr +  realinterest  +
                      +  stock_market +  old_age +  cit + 
                      +  std_net_lending_lag + emprot_reg_lag + factor( country), 
                    data = a.out$imputations$imp1) 

# Disply Table 1 
stargazer(container_basis, container_twoway, container_lag, 
          container_basis_emprot, container_twoway_emprot, container_lag_emprot,
          coef = list(mod1_ud_coef, mod2_ud_coef, mod3_ud_coef,
                      mod1_emprot_reg_coef,mod2_emprot_reg_coef, mod3_emprot_reg_coef),
          se = list(mod1_ud_se, mod2_ud_se, mod3_ud_se, 
                    mod1_emprot_reg_se, mod2_emprot_reg_se, mod3_emprot_reg_se), type = "text",
          covariate.labels = c("Trade Union Density",  "Employment Protection", "RTI Score", "FDI out (\\% GDP)", "Real GDP Growth", 
                               "Real Interests",  "Stock Capital.", "Old Age Dep.", "Corp. Income Tax",  "Net Lending Lag","Trade Union Density Lag",
                               "Employment Protecion Lag"),
          omit = c("year", "country", "Constant"), order=c(1,12,2,3,4,5,6,7,8,9,10),
          dep.var.labels = "Corporate Net Lending (\\% GDP)",  title = "Higher Trade Union Density and Stronger Employment Protection are associated with lower Corporate Savings",
          add.lines = list(c("Country Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                           c("Year Fixe Effects", "$\\times$", "$\\checkmark$", "$\\times$", "$\\times$", "$\\checkmark$", "$\\times$")))

#------------------------------------------------------------------------------------------------------------------------
# 4) Figure 2 - PARTILA RESIDUAL PLOTS 
#------------------------------------------------------------------------------------------------------------------------

model_ep<-  lm(std_net_lending ~ emprot_reg + factor(country) + factor(year) + realgdpgr + rti_score +
               FDI_standardized + realinterest +  old_age + stock_market + cit, data=a.out$imputations$imp1)

model_ud<-  lm(std_net_lending ~ ud + factor(country) +  factor(year) + realgdpgr + rti_score + 
                 FDI_standardized + realinterest + old_age + stock_market + cit, data=a.out$imputations$imp1)

visreg(model_ep, "emprot_reg" , gg=TRUE, 
       points.par = list(alpha = 0.7, size=1.5), line.par=list(col="grey25")) + theme_bw() + 
  theme(text=element_text(size=14, family="Lato-Light"))  + ylab("Corporate Savings (% GDP)") +
  xlab("Employment protection strictness") 

visreg(model_ud, "ud" , gg=TRUE, points.par = list(alpha = 0.7, size=1.5), line.par=list(col="grey25")) + theme_bw() + 
  theme(text=element_text(size=14, family="Lato-Light"))  + ylab("Corporate Savings (% GDP)") +
  xlab("Trade Union Density") 


#------------------------------------------------------------------------------------------------------------------------
# 4) APPENDIX
#------------------------------------------------------------------------------------------------------------------------


# FIGURE A.1: OVERIMPUTED VALUES IF TRADE UNION DENSITY 

plot(a.out, which.vars = "ud")
plot<-overimpute(a.out, var="ud", main = "", legend = T)

# TABLE A.2 SUMMARY STATISTICS 

summary<-subset(df, select=c("std_net_lending", "ud", "emprot_reg" ,"rti_score" , 
                             "FDI_standardized", "realgdpgr", "realinterest",
                             "stock_market", "old_age", "cit", "adjcov" ))

stargazer(summary, type="text", 
          covariate.labels = c("Corporate Net Lending (\\% GDP)",
                               "Trade Union Density", "Employment Protection", 
                               "RTI Score", "FDI out (% GDP)", "Real GDP Growth",
                               "Real Interests", "Stock Market Capitalization", "Old Age Dep.", 
                               "Corporate Income Tax", "Bargaining Coverage"))

# FIGURE A.2 DEVELOPMENT OF CORPORATE SAVINGS ACROSSS COUNTRIES 

ggplot(df, aes(year,std_net_lending)) + geom_line() + facet_wrap(~ country, scale="free")+ theme_bw() +
  theme(text=element_text(size=12, family="Frutiger-Light")) + xlab("") + ylab("Corporat Savings (% GDP)")

# FIGURE A.3 DEVELOPMENT OF TRADE UNION DENSITY ACROSSS COUNTRIES 

ggplot(df, aes(year, ud)) + geom_line() + facet_wrap(~country, scale="free")+ theme_bw() +
  theme(text=element_text(size=12, family="Frutiger-Light")) + xlab("") + ylab("Trade Union Density")

# FIGURE A.3 DEVELOPMENT OF EMPLOYMENT PROTECTION INDICATOR ACROSSS COUNTRIES 

ggplot(df, aes(year, emprot_reg)) + geom_line() + facet_wrap(~ country, scale="free")+ theme_bw() +
  theme(text=element_text(size=12, family="Frutiger-Light")) + xlab("") + ylab("Employment Protection Legislation")


# TABLE A.3 THE ASSOCIATION BETWEEN TRADE UNION DENSITY AND EMPLOYMENT PROTECTION LEGISLATION

for(i in 1:a.out$m) {}

# Table A.3 - Model 1

# Conduct regressions across all imputations 
mod1_ud <-lapply(a.out$imputations, function(x){
  panelAR(emprot_reg ~  ud + factor(country), 
          data = a.out$imputations[[i]],
          timeVar ="year",
          panelVar="country", autoCorr="psar1",
          panelCorrMethod="pcse")
})

# Extract standard errors
mod1_ud_se <- list()
mod1_ud_se <- lapply(1:5, function(x){
  sqrt(diag(mod1_ud[[x]]$vcov)) 
})

# Take averages of standard errors across the models 
mod1_ud_se <-as.numeric(Reduce("+", mod1_ud_se)/5)

# Extract coefficients
mod1_ud_coef <- list()
mod1_ud_coef <- lapply(1:5, function(x){
  mod1_ud[[x]]$coefficients
})

# Take averages of coefficients across the models 
mod1_ud_coef <- Reduce("+", mod1_ud_coef)/5


# Table A.3 - Model 2

# Conduct regressions across all imputations 
mod2_ud <-lapply(a.out$imputations, function(x){
  panelAR(emprot_reg ~  ud + rti_score +  FDI_standardized + realgdpgr + realinterest +  stock_market + old_age + cit + 
            + factor(country) + factor(year), 
          data = a.out$imputations[[i]],
          timeVar ="year",
          panelVar="country", autoCorr="psar1",
          panelCorrMethod="pcse")
})

# Extract standard errors
mod2_ud_se <- list()
mod2_ud_se <- lapply(1:5, function(x){
  sqrt(diag(mod2_ud[[x]]$vcov)) 
})

# Take averages of standard errors across the models 
mod2_ud_se <-as.numeric(Reduce("+", mod2_ud_se)/5)

# Extract coefficients
mod2_ud_coef <- list()
mod2_ud_coef <- lapply(1:5, function(x){
  mod2_ud[[x]]$coefficients
})

# Take averages of coefficients across the models 
mod2_ud_coef <- Reduce("+", mod2_ud_coef)/5

# Table A.3 - Model 3

mod3_emprot_reg <-lapply(a.out$imputations, function(x){
  pcse(lm(emprot_reg ~  ud  + rti_score + FDI_standardized  + realgdpgr + realinterest  +
            stock_market  + old_age + cit + std_net_lending_lag +  emprot_reg_lag + factor(country), data = a.out$imputations[[i]]),
       groupT =a.out$imputations[[i]]$year, groupN=a.out$imputations[[i]]$country)
})

# Extract standard errors
mod3_emprot_reg_se <- list()
mod3_emprot_reg_se <- lapply(1:5, function(x){
  mod3_emprot_reg[[x]]$pcse
})

# Take averages of standard errors across the models 
mod3_emprot_reg_se <-as.numeric(Reduce("+", mod3_emprot_reg_se)/5)

# Extract coefficients
mod3_emprot_reg_coef <- list()
mod3_emprot_reg_coef <- lapply(1:5, function(x){
  mod3_emprot_reg[[x]]$b
})

# Take averages of coefficients across the models 
mod3_emprot_reg_coef <- as.numeric( Reduce("+", mod3_emprot_reg_coef)/5)

## Prepare container models 

# Note: these are placeholder models to insert coeffiecents/ standard errors estimates above into an easy table
# with stargazer
container_basis_emprot<- lm( emprot_reg ~ ud + factor(country), 
                             data = a.out$imputations$imp1)  
container_twoway_emprot<- lm(emprot_reg ~  ud + rti_score +  FDI_standardized + realgdpgr + realinterest +  stock_market + old_age + cit + 
                               + factor(country) + factor(year), data = a.out$imputations$imp1)
container_lag_emprot<-lm(  emprot_reg ~  ud  +  rti_score +   FDI_standardized +
                             realgdpgr +  realinterest  +
                             +  stock_market +  old_age +  cit + 
                             +  ud_lag + emprot_reg_lag + factor( country), 
                           data = a.out$imputations$imp1) 

# Display Table A.3
stargazer(container_basis_emprot, container_twoway_emprot, container_lag_emprot,
          coef = list( mod1_ud_coef, mod2_ud_coef, mod3_emprot_reg_coef),
          se = list(mod1_ud_se, mod2_ud_se, mod3_emprot_reg_se), type = "text",
covariate.labels = c("Trade Union Density", "RTI Score", "FDI out (\\% GDP)", "Real GDP Growth", 
                               "Real Interests",  "Stock Capital.", "Old Age Dep.", "Corp. Income Tax",  "Employment Protection Lag","Trade Union Density Lag"),
          omit = c("year", "country", "Constant"), 
          dep.var.labels = "Employment Protection Legislation Index",  title = "Higher Trade Union Density is associated with Stronger Employment Protection Legislation",
          add.lines = list(c("Country Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                           c("Year-Fixe Effects", "$\\times$", "$\\checkmark$", "$\\times$")))


# TABLE A.4 TRADE UNION DENSITY AND CORPORATE SAVINGS - 5 YEAR PERIOD AVERAGES

# Average imputed data over 5 year period

imp_5y<-lapply(a.out$imputations, function(x){
  x <- x %>%
    arrange(country, year) %>%
    dplyr::group_by(country) %>%
    dplyr::mutate(period = gl(ceiling(n() / 5), 5, length = n())) %>%
    dplyr::group_by(country, period) %>%
    dplyr::summarise_all(mean) 
})


### Table A.4 - Model 1 

imp1 <- lapply(list(imp_5y$imp1, imp_5y$imp2, imp_5y$imp3, imp_5y$imp4, 
                    imp_5y$imp5),lm, formula=std_net_lending ~  ud + factor(country))

betas1 <- MIextract(imp1, fun = coef)
vars1 <- MIextract(imp1, fun = vcov)
MI_fixed_bi<-(summary(MIcombine(betas1, vars1)))
MI_fixed_bi_coef<-MI_fixed_bi$results
MI_fixed_bi_se<-MI_fixed_bi$se


### Table A.4 - Model 2

imp2 <- lapply(list(imp_5y$imp1, imp_5y$imp2,
                    imp_5y$imp3, imp_5y$imp4, 
                    imp_5y$imp2),lm, formula=std_net_lending ~ ud + rti_score + FDI_standardized + 
                 realgdpgr + realinterest +  stock_market + old_age + cit + 
                 + factor(country) + factor(period)) 

betas2 <- MIextract(imp2, fun = coef)
vars2 <- MIextract(imp2, fun = vcov)
MI_fixed_macro<-(summary(MIcombine(betas2, vars2)))
MI_fixed_macro_coef<-MI_fixed_macro$results
MI_fixed_macro_se<-MI_fixed_macro$se

### Table A.4 - Model 3

imp3 <- lapply(list(imp_5y$imp1, imp_5y$imp2,
                    imp_5y$imp3, imp_5y$imp4, 
                    imp_5y$imp2),lm, formula=std_net_lending ~  emprot_reg + factor(country))

betas3 <- MIextract(imp3, fun = coef)
vars3 <- MIextract(imp3, fun = vcov)
MI_fixed_bi_emprot<-(summary(MIcombine(betas3, vars3)))
MI_fixed_bi_emprot_coef<-MI_fixed_bi_emprot$results
MI_fixed_bi_emprot_se<-MI_fixed_bi$se


### Table A.4 - Model 4

imp4 <- lapply(list(imp_5y$imp1, imp_5y$imp2,
                    imp_5y$imp3, imp_5y$imp4, 
                    imp_5y$imp2),lm, formula=std_net_lending ~  emprot_reg + rti_score +  FDI_standardized + 
                 realgdpgr + realinterest +  stock_market + old_age + cit + 
                 + factor(country) )

betas4 <- MIextract(imp4, fun = coef)
vars4 <- MIextract(imp4, fun = vcov)
MI_fixed_macro<-(summary(MIcombine(betas4, vars4)))
MI_fixed_macro_emprot_coef<-MI_fixed_macro$results
MI_fixed_macro_emprot_se<-MI_fixed_macro$se


## Prepare container models 

# Note: these are placeholder models to insert coeffiecents/ standard errors estimates 
# above into an easy table with stargazer

container_basis<- lm( std_net_lending ~ ud + factor(country), 
                      data = imp_5y$imp1)  
container_twoway<- lm( std_net_lending ~  ud  +  rti_score +  FDI_standardized +  realgdpgr +
                         realinterest +    stock_market +  old_age +  cit +
                         factor( country) + factor( period), data = imp_5y$imp1)
container_basis_emprot<- lm( std_net_lending ~ emprot_reg + factor(country), 
                             data = imp_5y$imp1)  
container_twoway_emprot<- lm( std_net_lending ~  emprot_reg  +  rti_score +  FDI_standardized +  realgdpgr +
                                realinterest +    stock_market +  old_age +  cit +
                                factor(country) + factor(period), data = imp_5y$imp1)

stargazer(container_basis, container_twoway,
          container_basis_emprot, container_twoway_emprot, 
          coef = list(MI_fixed_bi_coef, MI_fixed_macro_coef, 
                      MI_fixed_bi_emprot_coef, MI_fixed_macro_emprot_coef),
          se = list(MI_fixed_bi_se, MI_fixed_macro_se, MI_fixed_bi_emprot_se, MI_fixed_macro_emprot_se ), type = "text",
          covariate.labels = c("Trade Union Density", "Employment Protection", "RTI Score", "FDI out (\\% GDP)", "Real GDP Growth", 
                               "Real Interests", "Stock Capital.", "Old Age Dep.", "Corp. Income Tax",  "Net Lending Lag"),
          omit = c("year", "country", "Constant", "period"), order=c(1,9,2,3,4,5,6,7,8),
          dep.var.labels = "Corporate Net Lending (\\% GDP)",  title = "Trade Union Density and Corporate Savings - 5 Year Period Averages",
          add.lines = list(c("Country Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                           c("Year Fixe Effects", "$\\times$", "$\\times$", "$\\times$", "$\\checkmark$", "$\\times$")))


# TABLE A.5 TRADE UNION DENSITY AND CORPORATE SAVINGS - 3 YEAR PERIOD AVERAGES

# Average imputed data over 3 year period

imp_3y<-lapply(a.out$imputations, function(x){
  x <- x %>%
    arrange(country, year) %>%
    dplyr::group_by(country) %>%
    dplyr::mutate(period = gl(ceiling(n() / 3), 3, length = n())) %>%
    dplyr::group_by(country, period) %>%
    dplyr::summarise_all(mean) 
})


### Table A.5 - Model 1 

imp1 <- lapply(list(imp_3y$imp1, imp_3y$imp2, imp_3y$imp3, imp_3y$imp4, 
                    imp_3y$imp5),lm, formula=std_net_lending ~  ud + factor(country))

betas1 <- MIextract(imp1, fun = coef)
vars1 <- MIextract(imp1, fun = vcov)
MI_fixed_bi<-(summary(MIcombine(betas1, vars1)))
MI_fixed_bi_coef<-MI_fixed_bi$results
MI_fixed_bi_se<-MI_fixed_bi$se


### Table A.5 - Model 2

imp2 <- lapply(list(imp_3y$imp1, imp_3y$imp2,
                    imp_3y$imp3, imp_3y$imp4, 
                    imp_3y$imp2),lm, formula=std_net_lending ~ ud + rti_score + FDI_standardized + 
                 realgdpgr + realinterest +  stock_market + old_age + cit + 
                 + factor(country) + factor(period)) 

betas2 <- MIextract(imp2, fun = coef)
vars2 <- MIextract(imp2, fun = vcov)
MI_fixed_macro<-(summary(MIcombine(betas2, vars2)))
MI_fixed_macro_coef<-MI_fixed_macro$results
MI_fixed_macro_se<-MI_fixed_macro$se

### Table A.5 - Model 3

imp3 <- lapply(list(imp_3y$imp1, imp_3y$imp2,
                    imp_3y$imp3, imp_3y$imp4, 
                    imp_3y$imp2),lm, formula=std_net_lending ~  emprot_reg + factor(country))

betas3 <- MIextract(imp3, fun = coef)
vars3 <- MIextract(imp3, fun = vcov)
MI_fixed_bi_emprot<-(summary(MIcombine(betas3, vars3)))
MI_fixed_bi_emprot_coef<-MI_fixed_bi_emprot$results
MI_fixed_bi_emprot_se<-MI_fixed_bi$se


### Table A.5 - Model 4

imp4 <- lapply(list(imp_3y$imp1, imp_3y$imp2,
                    imp_3y$imp3, imp_3y$imp4, 
                    imp_3y$imp2),lm, formula=std_net_lending ~  emprot_reg + rti_score +  FDI_standardized + 
                 realgdpgr + realinterest +  stock_market + old_age + cit + 
                 + factor(country) )

betas4 <- MIextract(imp4, fun = coef)
vars4 <- MIextract(imp4, fun = vcov)
MI_fixed_macro<-(summary(MIcombine(betas4, vars4)))
MI_fixed_macro_emprot_coef<-MI_fixed_macro$results
MI_fixed_macro_emprot_se<-MI_fixed_macro$se


## Prepare container models 

# Note: these are placeholder models to insert coeffiecents/ standard errors estimates 
# above into an easy table with stargazer

container_basis<- lm( std_net_lending ~ ud + factor(country), 
                      data = imp_3y$imp1)  
container_twoway<- lm( std_net_lending ~  ud  +  rti_score +  FDI_standardized +  realgdpgr +
                         realinterest +    stock_market +  old_age +  cit +
                         factor( country) + factor( period), data = imp_3y$imp1)
container_basis_emprot<- lm( std_net_lending ~ emprot_reg + factor(country), 
                             data = imp_5y$imp1)  
container_twoway_emprot<- lm( std_net_lending ~  emprot_reg  +  rti_score +  FDI_standardized +  realgdpgr +
                                realinterest +    stock_market +  old_age +  cit +
                                factor(country) + factor(period), data = imp_3y$imp1)

stargazer(container_basis, container_twoway,
          container_basis_emprot, container_twoway_emprot, 
          coef = list(MI_fixed_bi_coef, MI_fixed_macro_coef, 
                      MI_fixed_bi_emprot_coef, MI_fixed_macro_emprot_coef),
          se = list(MI_fixed_bi_se, MI_fixed_macro_se, MI_fixed_bi_emprot_se, MI_fixed_macro_emprot_se ), type = "text",
          covariate.labels = c("Trade Union Density", "Employment Protection", "RTI Score", "FDI out (\\% GDP)", "Real GDP Growth", 
                               "Real Interests", "Stock Capital.", "Old Age Dep.", "Corp. Income Tax",  "Net Lending Lag"),
          omit = c("year", "country", "Constant", "period"), order=c(1,9,2,3,4,5,6,7,8),
          dep.var.labels = "Corporate Net Lending (\\% GDP)",  title = "Trade Union Density and Corporate Savings - 3 Year Period Averages",
          add.lines = list(c("Country Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                           c("Year Fixe Effects", "$\\times$", "$\\times$", "$\\times$", "$\\checkmark$", "$\\times$")),
          out="/Users/nredekeer/Dropbox/02_PhD/07_Corporate_Savings/03Paper/10Resubmission2/tables/three_averages.tex")



# TABLE A.6 BARGAINING COVERAGE AND CORPORATE SAVINGS 

for(i in 1:a.out$m) {}

# Table A.6 - Model 1

# Conduct regressions across all imputations 
mod1_ud <-lapply(a.out$imputations, function(x){
  panelAR(std_net_lending ~  adjcov + factor(country), 
          data = a.out$imputations[[i]],
          timeVar ="year",
          panelVar="country", autoCorr="psar1",
          panelCorrMethod="pcse")
})

# Extract standard errors
mod1_ud_se <- list()
mod1_ud_se <- lapply(1:5, function(x){
  sqrt(diag(mod1_ud[[x]]$vcov)) 
})

# Take averages of standard errors across the models 
mod1_ud_se <-as.numeric(Reduce("+", mod1_ud_se)/5)

# Extract coefficients
mod1_ud_coef <- list()
mod1_ud_coef <- lapply(1:5, function(x){
  mod1_ud[[x]]$coefficients
})

# Take averages of coefficients across the models 
mod1_ud_coef <- Reduce("+", mod1_ud_coef)/5
df$crisis_dummy

# Table A.3 - Model 2

# Conduct regressions across all imputations 
mod2_ud <-lapply(a.out$imputations, function(x){
  panelAR(std_net_lending ~  adjcov + rti_score +  FDI_standardized + 
            realgdpgr + realinterest +  stock_market + old_age + cit + factor(country), 
          data = a.out$imputations[[i]],
          timeVar ="year",
          panelVar="country", autoCorr="psar1",
          panelCorrMethod="pcse")
})

# Extract standard errors
mod2_ud_se <- list()
mod2_ud_se <- lapply(1:5, function(x){
  sqrt(diag(mod2_ud[[x]]$vcov)) 
})

# Take averages of standard errors across the models 
mod2_ud_se <-as.numeric(Reduce("+", mod2_ud_se)/5)

# Extract coefficients
mod2_ud_coef <- list()
mod2_ud_coef <- lapply(1:5, function(x){
  mod2_ud[[x]]$coefficients
})

# Take averages of coefficients across the models 
mod2_ud_coef <- Reduce("+", mod2_ud_coef)/5

# Table A.3 - Model 3

mod3_emprot_reg <-lapply(a.out$imputations, function(x){
  pcse(lm(std_net_lending ~  adjcov  + rti_score + FDI_standardized  + realgdpgr + realinterest  +
            stock_market  + old_age + cit + factor(country) + factor(year), data = a.out$imputations[[i]]),
       groupT =a.out$imputations[[i]]$year, groupN=a.out$imputations[[i]]$country)
})

# Extract standard errors
mod3_emprot_reg_se <- list()
mod3_emprot_reg_se <- lapply(1:5, function(x){
  mod3_emprot_reg[[x]]$pcse
})

# Take averages of standard errors across the models 
mod3_emprot_reg_se <-as.numeric(Reduce("+", mod3_emprot_reg_se)/5)

# Extract coefficients
mod3_emprot_reg_coef <- list()
mod3_emprot_reg_coef <- lapply(1:5, function(x){
  mod3_emprot_reg[[x]]$b
})

# Take averages of coefficients across the models 
mod3_emprot_reg_coef <- as.numeric( Reduce("+", mod3_emprot_reg_coef)/5)

## Prepare container models 

# Note: these are placeholder models to insert coeffiecents/ standard errors estimates above into an easy table
# with stargazer
container_basis_emprot<- lm( std_net_lending ~ adjcov + factor(country), 
                             data = a.out$imputations$imp1)  
container_twoway_emprot<- lm(std_net_lending ~  adjcov + rti_score +  FDI_standardized + realgdpgr + realinterest +  stock_market + old_age + cit + 
                               + factor(country), data = a.out$imputations$imp1)
container_lag_emprot<-lm(  std_net_lending ~  adjcov  +  rti_score +   FDI_standardized +
                             realgdpgr +  realinterest  +
                             +  stock_market +  old_age +  cit + factor( country) + factor(year), 
                           data = a.out$imputations$imp1) 

# Display Table A.3
stargazer(container_basis_emprot, container_twoway_emprot, container_lag_emprot,
          coef = list( mod1_ud_coef, mod2_ud_coef, mod3_emprot_reg_coef),
          se = list(mod1_ud_se, mod2_ud_se, mod3_emprot_reg_se), type = "text",
          covariate.labels = c("Bargaining Coverage", "RTI Score", "FDI out (\\% GDP)", "Real GDP Growth", 
                               "Real Interests",  "Stock Capital.", "Old Age Dep.", "Corp. Income Tax"),
          omit = c("year", "country", "Constant"), 
          dep.var.labels = "Corporate Savings (% of GDP)", 
          title = "Higher Bargaining Coverage is associated with lower Corporate Savings",
          add.lines = list(c("Country Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                           c("Year-Fixe Effects", "$\\times$", "$\\times$", "$\\checkmark$")),
          out="/Users/nredekeer/Dropbox/02_PhD/07_Corporate_Savings/03Paper/10Resubmission2/tables/cross_counrty_alt.tex")






